{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Batch Correction for TCGA DNA Methylation Data\n",
    "In this notebook, I want to remove batch effects from the previously downloaded DNA methylation data.\n",
    "The obtained data was already preprocessed such that we have two sample matrices $S_{tumor} \\in \\mathbb{R}^{N \\times M_1}$ and $S_{normal} \\in \\mathbb{R}^{N \\times M_2}$ with $M_1$ and $M_2$ being the number of samples for tumor and normal tissues across all cancer types and $N$ being the number of genes.\n",
    "\n",
    "**That is, we already have computed the average promoter DNA methylation across all measured CpG sites per gene.**\n",
    "\n",
    "Next, we want to remove the batch effects using *ComBat*. As batches, we use the plate IDs as suggested in multiple articles.\n",
    "\n",
    "The workflow of that notebook is as follows:\n",
    "1. Load the big sample matrices for tumor and normal\n",
    "2. Split them into smaller gene-sample matrices for each cancer type and write them to disk\n",
    "3. Call a R script which does the batch correction with ComBat\n",
    "\n",
    "The results from batch correction are then read by another script and the final feature matrix is computed there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/sasse/.local/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "from matplotlib import gridspec\n",
    "plt.rc('font', family='Times New Roman')\n",
    "import h5py\n",
    "import seaborn as sns\n",
    "import os\n",
    "\n",
    "from sklearn.manifold import TSNE\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load Tumor and Normal Sample Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TCGA-F2-6880-01A-11D-2157-05|paad|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-2E-A9G8-01A-11D-A409-05|ucec|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-DQ-7588-01A-11D-2079-05|hnsc|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-78-7166-01A-12D-2064-05|luad|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-CQ-7068-01A-11D-2079-05|hnsc|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-P3-A6T4-01A-11D-A34K-05|hnsc|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-CX-7082-01A-11D-2014-05|hnsc|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-55-8621-01A-11D-2398-05|luad|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-BA-4077-01B-01D-1433-05|hnsc|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-44-7671-01A-11D-2064-05|luad|mean_beta_value_promoter</th>\n",
       "      <th>...</th>\n",
       "      <th>TCGA-J4-AATZ-01A-11D-A41L-05|prad|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-ET-A39N-01A-11D-A19K-05|thca|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-5L-AAT0-01A-12D-A41Q-05|brca|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-G2-A3VY-01A-11D-A231-05|blca|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-G9-6338-01A-12D-1963-05|prad|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-BS-A0V4-01A-11D-A14H-05|ucec|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-B0-4819-01A-01D-1275-05|kirc|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-06-0210-01A-01D-A45W-05|gbm|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-44-2665-01A-01D-A276-05|luad|mean_beta_value_promoter</th>\n",
       "      <th>TCGA-D1-A17K-01A-11D-A12K-05|ucec|mean_beta_value_promoter</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Symbol</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A1BG</th>\n",
       "      <td>0.443866</td>\n",
       "      <td>0.320773</td>\n",
       "      <td>0.723260</td>\n",
       "      <td>0.531928</td>\n",
       "      <td>0.908139</td>\n",
       "      <td>0.780224</td>\n",
       "      <td>0.405876</td>\n",
       "      <td>0.717752</td>\n",
       "      <td>0.736370</td>\n",
       "      <td>0.600592</td>\n",
       "      <td>...</td>\n",
       "      <td>0.557431</td>\n",
       "      <td>0.802080</td>\n",
       "      <td>0.680144</td>\n",
       "      <td>0.212928</td>\n",
       "      <td>0.735104</td>\n",
       "      <td>0.467060</td>\n",
       "      <td>0.709056</td>\n",
       "      <td>0.739527</td>\n",
       "      <td>0.628191</td>\n",
       "      <td>0.693468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A1BG-AS1</th>\n",
       "      <td>0.110391</td>\n",
       "      <td>0.025817</td>\n",
       "      <td>0.050617</td>\n",
       "      <td>0.073089</td>\n",
       "      <td>0.116924</td>\n",
       "      <td>0.036841</td>\n",
       "      <td>0.216182</td>\n",
       "      <td>0.100127</td>\n",
       "      <td>0.104327</td>\n",
       "      <td>0.077859</td>\n",
       "      <td>...</td>\n",
       "      <td>0.091708</td>\n",
       "      <td>0.049971</td>\n",
       "      <td>0.041602</td>\n",
       "      <td>0.084231</td>\n",
       "      <td>0.197002</td>\n",
       "      <td>0.211159</td>\n",
       "      <td>0.098477</td>\n",
       "      <td>0.033417</td>\n",
       "      <td>0.059654</td>\n",
       "      <td>0.045232</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A1CF</th>\n",
       "      <td>0.391441</td>\n",
       "      <td>0.250301</td>\n",
       "      <td>0.462136</td>\n",
       "      <td>0.567162</td>\n",
       "      <td>0.695086</td>\n",
       "      <td>0.543520</td>\n",
       "      <td>0.409346</td>\n",
       "      <td>0.785829</td>\n",
       "      <td>0.622589</td>\n",
       "      <td>0.663200</td>\n",
       "      <td>...</td>\n",
       "      <td>0.359442</td>\n",
       "      <td>0.891064</td>\n",
       "      <td>0.675270</td>\n",
       "      <td>0.363130</td>\n",
       "      <td>0.436892</td>\n",
       "      <td>0.403012</td>\n",
       "      <td>0.715297</td>\n",
       "      <td>0.700912</td>\n",
       "      <td>0.734370</td>\n",
       "      <td>0.234995</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A2M</th>\n",
       "      <td>0.810426</td>\n",
       "      <td>0.464732</td>\n",
       "      <td>0.529455</td>\n",
       "      <td>0.606023</td>\n",
       "      <td>0.491957</td>\n",
       "      <td>0.595080</td>\n",
       "      <td>0.445457</td>\n",
       "      <td>0.499083</td>\n",
       "      <td>0.557930</td>\n",
       "      <td>0.581326</td>\n",
       "      <td>...</td>\n",
       "      <td>0.295140</td>\n",
       "      <td>0.534252</td>\n",
       "      <td>0.463106</td>\n",
       "      <td>0.732270</td>\n",
       "      <td>0.574303</td>\n",
       "      <td>0.522476</td>\n",
       "      <td>0.534992</td>\n",
       "      <td>0.205325</td>\n",
       "      <td>0.434673</td>\n",
       "      <td>0.445305</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A2ML1</th>\n",
       "      <td>0.703449</td>\n",
       "      <td>0.297757</td>\n",
       "      <td>0.311103</td>\n",
       "      <td>0.529984</td>\n",
       "      <td>0.415163</td>\n",
       "      <td>0.499372</td>\n",
       "      <td>0.320599</td>\n",
       "      <td>0.727881</td>\n",
       "      <td>0.397720</td>\n",
       "      <td>0.674001</td>\n",
       "      <td>...</td>\n",
       "      <td>0.332043</td>\n",
       "      <td>0.778914</td>\n",
       "      <td>0.792049</td>\n",
       "      <td>0.192367</td>\n",
       "      <td>0.347380</td>\n",
       "      <td>0.310234</td>\n",
       "      <td>0.803575</td>\n",
       "      <td>0.751523</td>\n",
       "      <td>0.708201</td>\n",
       "      <td>0.415714</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 7391 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          TCGA-F2-6880-01A-11D-2157-05|paad|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.443866            \n",
       "A1BG-AS1                                           0.110391            \n",
       "A1CF                                               0.391441            \n",
       "A2M                                                0.810426            \n",
       "A2ML1                                              0.703449            \n",
       "\n",
       "          TCGA-2E-A9G8-01A-11D-A409-05|ucec|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.320773            \n",
       "A1BG-AS1                                           0.025817            \n",
       "A1CF                                               0.250301            \n",
       "A2M                                                0.464732            \n",
       "A2ML1                                              0.297757            \n",
       "\n",
       "          TCGA-DQ-7588-01A-11D-2079-05|hnsc|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.723260            \n",
       "A1BG-AS1                                           0.050617            \n",
       "A1CF                                               0.462136            \n",
       "A2M                                                0.529455            \n",
       "A2ML1                                              0.311103            \n",
       "\n",
       "          TCGA-78-7166-01A-12D-2064-05|luad|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.531928            \n",
       "A1BG-AS1                                           0.073089            \n",
       "A1CF                                               0.567162            \n",
       "A2M                                                0.606023            \n",
       "A2ML1                                              0.529984            \n",
       "\n",
       "          TCGA-CQ-7068-01A-11D-2079-05|hnsc|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.908139            \n",
       "A1BG-AS1                                           0.116924            \n",
       "A1CF                                               0.695086            \n",
       "A2M                                                0.491957            \n",
       "A2ML1                                              0.415163            \n",
       "\n",
       "          TCGA-P3-A6T4-01A-11D-A34K-05|hnsc|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.780224            \n",
       "A1BG-AS1                                           0.036841            \n",
       "A1CF                                               0.543520            \n",
       "A2M                                                0.595080            \n",
       "A2ML1                                              0.499372            \n",
       "\n",
       "          TCGA-CX-7082-01A-11D-2014-05|hnsc|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.405876            \n",
       "A1BG-AS1                                           0.216182            \n",
       "A1CF                                               0.409346            \n",
       "A2M                                                0.445457            \n",
       "A2ML1                                              0.320599            \n",
       "\n",
       "          TCGA-55-8621-01A-11D-2398-05|luad|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.717752            \n",
       "A1BG-AS1                                           0.100127            \n",
       "A1CF                                               0.785829            \n",
       "A2M                                                0.499083            \n",
       "A2ML1                                              0.727881            \n",
       "\n",
       "          TCGA-BA-4077-01B-01D-1433-05|hnsc|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.736370            \n",
       "A1BG-AS1                                           0.104327            \n",
       "A1CF                                               0.622589            \n",
       "A2M                                                0.557930            \n",
       "A2ML1                                              0.397720            \n",
       "\n",
       "          TCGA-44-7671-01A-11D-2064-05|luad|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.600592            \n",
       "A1BG-AS1                                           0.077859            \n",
       "A1CF                                               0.663200            \n",
       "A2M                                                0.581326            \n",
       "A2ML1                                              0.674001            \n",
       "\n",
       "                                     ...                              \\\n",
       "Symbol                               ...                               \n",
       "A1BG                                 ...                               \n",
       "A1BG-AS1                             ...                               \n",
       "A1CF                                 ...                               \n",
       "A2M                                  ...                               \n",
       "A2ML1                                ...                               \n",
       "\n",
       "          TCGA-J4-AATZ-01A-11D-A41L-05|prad|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.557431            \n",
       "A1BG-AS1                                           0.091708            \n",
       "A1CF                                               0.359442            \n",
       "A2M                                                0.295140            \n",
       "A2ML1                                              0.332043            \n",
       "\n",
       "          TCGA-ET-A39N-01A-11D-A19K-05|thca|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.802080            \n",
       "A1BG-AS1                                           0.049971            \n",
       "A1CF                                               0.891064            \n",
       "A2M                                                0.534252            \n",
       "A2ML1                                              0.778914            \n",
       "\n",
       "          TCGA-5L-AAT0-01A-12D-A41Q-05|brca|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.680144            \n",
       "A1BG-AS1                                           0.041602            \n",
       "A1CF                                               0.675270            \n",
       "A2M                                                0.463106            \n",
       "A2ML1                                              0.792049            \n",
       "\n",
       "          TCGA-G2-A3VY-01A-11D-A231-05|blca|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.212928            \n",
       "A1BG-AS1                                           0.084231            \n",
       "A1CF                                               0.363130            \n",
       "A2M                                                0.732270            \n",
       "A2ML1                                              0.192367            \n",
       "\n",
       "          TCGA-G9-6338-01A-12D-1963-05|prad|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.735104            \n",
       "A1BG-AS1                                           0.197002            \n",
       "A1CF                                               0.436892            \n",
       "A2M                                                0.574303            \n",
       "A2ML1                                              0.347380            \n",
       "\n",
       "          TCGA-BS-A0V4-01A-11D-A14H-05|ucec|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.467060            \n",
       "A1BG-AS1                                           0.211159            \n",
       "A1CF                                               0.403012            \n",
       "A2M                                                0.522476            \n",
       "A2ML1                                              0.310234            \n",
       "\n",
       "          TCGA-B0-4819-01A-01D-1275-05|kirc|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.709056            \n",
       "A1BG-AS1                                           0.098477            \n",
       "A1CF                                               0.715297            \n",
       "A2M                                                0.534992            \n",
       "A2ML1                                              0.803575            \n",
       "\n",
       "          TCGA-06-0210-01A-01D-A45W-05|gbm|mean_beta_value_promoter  \\\n",
       "Symbol                                                                \n",
       "A1BG                                               0.739527           \n",
       "A1BG-AS1                                           0.033417           \n",
       "A1CF                                               0.700912           \n",
       "A2M                                                0.205325           \n",
       "A2ML1                                              0.751523           \n",
       "\n",
       "          TCGA-44-2665-01A-01D-A276-05|luad|mean_beta_value_promoter  \\\n",
       "Symbol                                                                 \n",
       "A1BG                                               0.628191            \n",
       "A1BG-AS1                                           0.059654            \n",
       "A1CF                                               0.734370            \n",
       "A2M                                                0.434673            \n",
       "A2ML1                                              0.708201            \n",
       "\n",
       "          TCGA-D1-A17K-01A-11D-A12K-05|ucec|mean_beta_value_promoter  \n",
       "Symbol                                                                \n",
       "A1BG                                               0.693468           \n",
       "A1BG-AS1                                           0.045232           \n",
       "A1CF                                               0.234995           \n",
       "A2M                                                0.445305           \n",
       "A2ML1                                              0.415714           \n",
       "\n",
       "[5 rows x 7391 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tumor_samples = pd.read_csv('../../data/pancancer/TCGA/methylation/download_new/tumor_sample_matrix.tsv.gz',\n",
    "                               compression='gzip',\n",
    "                               sep='\\t')\n",
    "# rename first column to 'Symbol'\n",
    "tumor_samples.columns = ['Symbol'] + tumor_samples.columns[1:].tolist()\n",
    "# get rid of all the support columns\n",
    "tumor_betaval_cols = ['Symbol'] + [i for i in tumor_samples.columns[1:] if \"mean_beta_value_promoter\" in i]\n",
    "tumor_samples = tumor_samples[tumor_betaval_cols]\n",
    "# put the symbol as index\n",
    "tumor_samples.set_index('Symbol', inplace=True)\n",
    "# get cancer type, transform and add it as column\n",
    "ctype_tumor = [i.split('|')[1].upper() for i in tumor_samples.columns]\n",
    "tumor_samples_t = tumor_samples.T\n",
    "tumor_samples_t['ctype'] = ctype_tumor\n",
    "tumor_samples.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Symbol</th>\n",
       "      <th>A1BG</th>\n",
       "      <th>A1BG-AS1</th>\n",
       "      <th>A1CF</th>\n",
       "      <th>A2M</th>\n",
       "      <th>A2ML1</th>\n",
       "      <th>A4GNT</th>\n",
       "      <th>AAAS</th>\n",
       "      <th>AACS</th>\n",
       "      <th>AADAC</th>\n",
       "      <th>AADACL2</th>\n",
       "      <th>...</th>\n",
       "      <th>ZXDA</th>\n",
       "      <th>ZXDB</th>\n",
       "      <th>ZXDC</th>\n",
       "      <th>ZYG11A</th>\n",
       "      <th>ZYG11B</th>\n",
       "      <th>ZYX</th>\n",
       "      <th>ZZEF1</th>\n",
       "      <th>ZZZ3</th>\n",
       "      <th>yR211F11.2</th>\n",
       "      <th>ctype</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_value_promoter</th>\n",
       "      <td>0.624236</td>\n",
       "      <td>0.044528</td>\n",
       "      <td>0.838009</td>\n",
       "      <td>0.598340</td>\n",
       "      <td>0.761133</td>\n",
       "      <td>0.693980</td>\n",
       "      <td>0.178805</td>\n",
       "      <td>0.063790</td>\n",
       "      <td>0.784948</td>\n",
       "      <td>0.898372</td>\n",
       "      <td>...</td>\n",
       "      <td>0.068541</td>\n",
       "      <td>0.187936</td>\n",
       "      <td>0.076922</td>\n",
       "      <td>0.291098</td>\n",
       "      <td>0.060295</td>\n",
       "      <td>0.047897</td>\n",
       "      <td>0.075362</td>\n",
       "      <td>0.033015</td>\n",
       "      <td>0.919028</td>\n",
       "      <td>LUSC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_value_promoter</th>\n",
       "      <td>0.796221</td>\n",
       "      <td>0.060968</td>\n",
       "      <td>0.853203</td>\n",
       "      <td>0.426811</td>\n",
       "      <td>0.751024</td>\n",
       "      <td>0.773104</td>\n",
       "      <td>0.205047</td>\n",
       "      <td>0.068452</td>\n",
       "      <td>0.889046</td>\n",
       "      <td>0.889440</td>\n",
       "      <td>...</td>\n",
       "      <td>0.490259</td>\n",
       "      <td>0.453815</td>\n",
       "      <td>0.110322</td>\n",
       "      <td>0.338830</td>\n",
       "      <td>0.072566</td>\n",
       "      <td>0.056229</td>\n",
       "      <td>0.086861</td>\n",
       "      <td>0.036836</td>\n",
       "      <td>0.925634</td>\n",
       "      <td>BLCA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_value_promoter</th>\n",
       "      <td>0.703644</td>\n",
       "      <td>0.061137</td>\n",
       "      <td>0.627578</td>\n",
       "      <td>0.530384</td>\n",
       "      <td>0.613380</td>\n",
       "      <td>0.551860</td>\n",
       "      <td>0.234964</td>\n",
       "      <td>0.047511</td>\n",
       "      <td>0.865786</td>\n",
       "      <td>0.905003</td>\n",
       "      <td>...</td>\n",
       "      <td>0.130354</td>\n",
       "      <td>0.154290</td>\n",
       "      <td>0.082387</td>\n",
       "      <td>0.404636</td>\n",
       "      <td>0.054379</td>\n",
       "      <td>0.041391</td>\n",
       "      <td>0.070617</td>\n",
       "      <td>0.032585</td>\n",
       "      <td>0.922868</td>\n",
       "      <td>COAD</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_value_promoter</th>\n",
       "      <td>0.745175</td>\n",
       "      <td>0.051093</td>\n",
       "      <td>0.810508</td>\n",
       "      <td>0.222311</td>\n",
       "      <td>0.834544</td>\n",
       "      <td>0.742291</td>\n",
       "      <td>0.181212</td>\n",
       "      <td>0.056555</td>\n",
       "      <td>0.837956</td>\n",
       "      <td>0.837963</td>\n",
       "      <td>...</td>\n",
       "      <td>0.095871</td>\n",
       "      <td>0.245532</td>\n",
       "      <td>0.088191</td>\n",
       "      <td>0.304986</td>\n",
       "      <td>0.061216</td>\n",
       "      <td>0.056205</td>\n",
       "      <td>0.083569</td>\n",
       "      <td>0.043757</td>\n",
       "      <td>0.882693</td>\n",
       "      <td>BLCA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_value_promoter</th>\n",
       "      <td>0.684194</td>\n",
       "      <td>0.065420</td>\n",
       "      <td>0.662524</td>\n",
       "      <td>0.522890</td>\n",
       "      <td>0.730755</td>\n",
       "      <td>0.685777</td>\n",
       "      <td>0.204331</td>\n",
       "      <td>0.091474</td>\n",
       "      <td>0.683887</td>\n",
       "      <td>0.823476</td>\n",
       "      <td>...</td>\n",
       "      <td>0.474710</td>\n",
       "      <td>0.450910</td>\n",
       "      <td>0.117562</td>\n",
       "      <td>0.350380</td>\n",
       "      <td>0.085787</td>\n",
       "      <td>0.063076</td>\n",
       "      <td>0.074156</td>\n",
       "      <td>0.067302</td>\n",
       "      <td>0.876185</td>\n",
       "      <td>LUAD</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 28353 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "Symbol                                                  A1BG  A1BG-AS1  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.624236  0.044528   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.796221  0.060968   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.703644  0.061137   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.745175  0.051093   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.684194  0.065420   \n",
       "\n",
       "Symbol                                                  A1CF       A2M  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.838009  0.598340   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.853203  0.426811   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.627578  0.530384   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.810508  0.222311   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.662524  0.522890   \n",
       "\n",
       "Symbol                                                 A2ML1     A4GNT  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.761133  0.693980   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.751024  0.773104   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.613380  0.551860   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.834544  0.742291   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.730755  0.685777   \n",
       "\n",
       "Symbol                                                  AAAS      AACS  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.178805  0.063790   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.205047  0.068452   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.234964  0.047511   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.181212  0.056555   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.204331  0.091474   \n",
       "\n",
       "Symbol                                                 AADAC   AADACL2  ...    \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.784948  0.898372  ...     \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.889046  0.889440  ...     \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.865786  0.905003  ...     \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.837956  0.837963  ...     \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.683887  0.823476  ...     \n",
       "\n",
       "Symbol                                                  ZXDA      ZXDB  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.068541  0.187936   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.490259  0.453815   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.130354  0.154290   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.095871  0.245532   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.474710  0.450910   \n",
       "\n",
       "Symbol                                                  ZXDC    ZYG11A  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.076922  0.291098   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.110322  0.338830   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.082387  0.404636   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.088191  0.304986   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.117562  0.350380   \n",
       "\n",
       "Symbol                                                ZYG11B       ZYX  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.060295  0.047897   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.072566  0.056229   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.054379  0.041391   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.061216  0.056205   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.085787  0.063076   \n",
       "\n",
       "Symbol                                                 ZZEF1      ZZZ3  \\\n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...  0.075362  0.033015   \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...  0.086861  0.036836   \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...  0.070617  0.032585   \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...  0.083569  0.043757   \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...  0.074156  0.067302   \n",
       "\n",
       "Symbol                                              yR211F11.2  ctype  \n",
       "TCGA-22-5471-11A-01D-1633-05|lusc|mean_beta_val...    0.919028   LUSC  \n",
       "TCGA-GC-A3WC-11A-11D-A231-05|blca|mean_beta_val...    0.925634   BLCA  \n",
       "TCGA-AA-3713-11A-01D-1721-05|coad|mean_beta_val...    0.922868   COAD  \n",
       "TCGA-CU-A0YN-11A-11D-A10W-05|blca|mean_beta_val...    0.882693   BLCA  \n",
       "TCGA-38-4631-11A-01D-1756-05|luad|mean_beta_val...    0.876185   LUAD  \n",
       "\n",
       "[5 rows x 28353 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normal_samples = pd.read_csv('../../data/pancancer/TCGA/methylation/download_new/normal_sample_matrix.tsv.gz',\n",
    "                             compression='gzip',\n",
    "                             sep='\\t')\n",
    "# put symbol in first column name\n",
    "normal_samples.columns = ['Symbol'] + normal_samples.columns[1:].tolist()\n",
    "# get rid of all the support columns\n",
    "normal_betaval_cols = ['Symbol'] + [i for i in normal_samples.columns[1:] if \"mean_beta_value_promoter\" in i]\n",
    "normal_samples = normal_samples[normal_betaval_cols]\n",
    "# put the symbol as index\n",
    "normal_samples.set_index('Symbol', inplace=True)\n",
    "# get cancer type, transform and add it as column\n",
    "ctype_normal = [i.split('|')[1].upper() for i in normal_samples.columns]\n",
    "normal_samples_t = normal_samples.T\n",
    "normal_samples_t['ctype'] = ctype_normal\n",
    "normal_samples_t.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Split per Cancer Type and write small matrices to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(28352, 419)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cols_of_interest = [i for i in tumor_samples.columns[1:] if i.split('|')[1] == 'blca']\n",
    "tumor_samples[cols_of_interest].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote matrices for cancer type PAAD\n",
      "Wrote matrices for cancer type UCEC\n",
      "Wrote matrices for cancer type HNSC\n",
      "Wrote matrices for cancer type LUAD\n",
      "Wrote matrices for cancer type KIRP\n",
      "Wrote matrices for cancer type KIRC\n",
      "Wrote matrices for cancer type PRAD\n",
      "Wrote matrices for cancer type COAD\n",
      "Wrote matrices for cancer type GBM\n",
      "Wrote matrices for cancer type LUSC\n",
      "Wrote matrices for cancer type CESC\n",
      "Wrote matrices for cancer type SARC\n",
      "Wrote matrices for cancer type STAD\n",
      "Wrote matrices for cancer type BRCA\n",
      "Wrote matrices for cancer type BLCA\n",
      "Wrote matrices for cancer type THCA\n",
      "Wrote matrices for cancer type LIHC\n",
      "Wrote matrices for cancer type SKCM\n",
      "Wrote matrices for cancer type READ\n",
      "Wrote matrices for cancer type ESCA\n",
      "Wrote matrices for cancer type OV\n"
     ]
    }
   ],
   "source": [
    "base_dir = '../../data/pancancer/TCGA/methylation/gene_sample_matrices_cancertype'\n",
    "\n",
    "processed_cancertypes = []\n",
    "for ctype in tumor_samples_t.ctype.unique():\n",
    "    # tumor samples\n",
    "    cols_of_interest = [i for i in tumor_samples.columns if i.split('|')[1].upper() == ctype]\n",
    "    tumor_samples_ctype = tumor_samples[cols_of_interest]\n",
    "    \n",
    "    # normal samples\n",
    "    cols_of_interest = [i for i in normal_samples.columns if i.split('|')[1].upper() == ctype]\n",
    "    normal_samples_ctype = normal_samples[cols_of_interest]\n",
    "    \n",
    "    # don't process the cancer type when there are no samples\n",
    "    if tumor_samples_ctype.shape[1] <= 1 or normal_samples_ctype.shape[1] <= 1:\n",
    "        print (\"Warning: No normal or tumor data for {}... Not processing it\".format(ctype))\n",
    "        continue\n",
    "\n",
    "    # directory\n",
    "    ctype_dir = os.path.join(base_dir, ctype)\n",
    "    if not os.path.isdir(ctype_dir):\n",
    "        os.mkdir(ctype_dir)\n",
    "\n",
    "    # write them to disk\n",
    "    tumor_samples_ctype.dropna(axis=0).to_csv(os.path.join(ctype_dir, 'tumor_samples.tsv'), sep='\\t')\n",
    "    normal_samples_ctype.dropna(axis=0).to_csv(os.path.join(ctype_dir, 'normal_samples.tsv'), sep='\\t')\n",
    "    \n",
    "    # pheno data\n",
    "    pheno_t = pd.DataFrame(tumor_samples_ctype.columns, columns=['Name'])\n",
    "    pheno_t['index'] = np.arange(pheno_t.shape[0])\n",
    "    pheno_t['cancer'] = 'tumor'\n",
    "    pheno_t['batch'] = [i[0].split('-')[5] for i in pheno_t.Name.str.split('|')]\n",
    "    pheno_t.set_index('Name', inplace=True)\n",
    "    pheno_t.to_csv(os.path.join(ctype_dir, 'pheno_tumor.tsv'), sep='\\t')\n",
    "\n",
    "    pheno_n = pd.DataFrame(normal_samples_ctype.columns, columns=['Name'])\n",
    "    pheno_n['index'] = np.arange(pheno_n.shape[0])\n",
    "    pheno_n['cancer'] = 'normal'\n",
    "    pheno_n['batch'] = [i[0].split('-')[5] for i in pheno_n.Name.str.split('|')]\n",
    "    pheno_n.set_index('Name', inplace=True)\n",
    "    pheno_n.to_csv(os.path.join(ctype_dir, 'pheno_normal.tsv'), sep='\\t')\n",
    "\n",
    "    print (\"Wrote matrices for cancer type {}\".format(ctype))\n",
    "    processed_cancertypes.append(ctype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Do batch correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processed PAAD\n",
      "Processed UCEC\n",
      "Processed HNSC\n",
      "Processed LUAD\n",
      "Processed KIRP\n",
      "Processed KIRC\n",
      "Processed PRAD\n",
      "Processed COAD\n",
      "Processed GBM\n",
      "Processed LUSC\n",
      "Processed CESC\n",
      "Processed SARC\n",
      "Processed STAD\n",
      "Processed BRCA\n",
      "Processed BLCA\n",
      "Processed THCA\n",
      "Processed LIHC\n",
      "Processed SKCM\n",
      "Processed READ\n",
      "Processed ESCA\n"
     ]
    }
   ],
   "source": [
    "import subprocess\n",
    "\n",
    "base_path = '../../data/pancancer/TCGA/methylation/gene_sample_matrices_cancertype/{}'\n",
    "call = 'Rscript batch_correction.R {} {} {}'\n",
    "for ctype in processed_cancertypes:\n",
    "    ctype_dir = base_path.format(ctype)\n",
    "    tumor_sample_path = os.path.join(ctype_dir, 'tumor_samples.tsv')\n",
    "    tumor_pheno_path = os.path.join(ctype_dir, 'pheno_tumor.tsv')\n",
    "    tumor_out_path = os.path.join(ctype_dir, 'tumor_samples.adjusted.tsv')\n",
    "    subprocess.call(call.format(tumor_sample_path, tumor_pheno_path, tumor_out_path), shell=True)\n",
    "\n",
    "    normal_sample_path = os.path.join(ctype_dir, 'normal_samples.tsv')\n",
    "    normal_pheno_path = os.path.join(ctype_dir, 'pheno_normal.tsv')\n",
    "    normal_out_path = os.path.join(ctype_dir, 'normal_samples.adjusted.tsv')\n",
    "    subprocess.call(call.format(normal_sample_path, normal_pheno_path, normal_out_path), shell=True)\n",
    "    print (\"Processed {}\".format(ctype))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Plot a PCA to see how well the batches cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 2)"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "from sklearn.manifold import TSNE\n",
    "emb = TSNE(n_components=2).fit_transform(tumor_samples_ctype.dropna(axis=0).T)\n",
    "emb.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(20, 10))\n",
    "cmap = pd.DataFrame(pd.Series([i[0].split('-')[5] for i in tumor_samples_ctype.columns[1:].str.split('|')]).unique(), columns=['ID'])\n",
    "cmap['number'] = np.arange(cmap.shape[0])\n",
    "cmap.set_index('ID', inplace=True)\n",
    "codes = pd.DataFrame([i[0].split('-')[5] for i in tumor_samples_ctype.columns.str.split('|')], columns=['ID'])\n",
    "colors = codes.join(mapping, on='ID').number\n",
    "plt.scatter(emb[:, 0], emb[:, 1], c=colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'colors' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-62-987f35b8a039>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mfig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0memb_new\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTSNE\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_components\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0madjusted_tumor\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0memb_new\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0memb_new\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcolors\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'colors' is not defined"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1440x720 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "adjusted_tumor = pd.read_csv('../../data/pancancer/TCGA/methylation/gene_sample_matrices_cancertype/PAAD/tumor_samples.adjusted.tsv', sep='\\t')\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "emb_new = TSNE(n_components=2).fit_transform(adjusted_tumor.T)\n",
    "plt.scatter(emb_new[:, 0], emb_new[:, 1], c=colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
